CSM - Lab4 - Prof. Minami

Codificação de Imagem com DWT

A. Objetivos:

  1. Efetuar a Codificação de Imagem e a Decodificação por DWT e IDWT
  2. Testar funções de Codificação Multinível
  3. Verificar a taxa de compressão só com a Componente de Aproximação
In [66]:
import numpy as np
import cv2 as cv
import pywt
import pywt.data
import matplotlib.image as mpimg
import matplotlib.pyplot as plt

B. Monte o seu google drive e verifique os seus arquivos:

In [67]:
from google.colab import drive
drive.mount('/content/gdrive')

Drive already mounted at /content/gdrive; to attempt to forcibly remount, call drive.mount("/content/gdrive", force_remount=True).

C. Codificação de Luminância (P&B) com DWT para a Pimentas

In [68]:
#img = mpimg.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/lucas.jpeg')
img = mpimg.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno.png')
# Conversão de BGRA (formato original da imagem) para RGB
img = cv.cvtColor(img, cv.COLOR_BGRA2BGR)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)


coefs2 = pywt.dwt2(img_gray,'haar', mode='periodization')  #1 nível de DWT
(cA, (cH, cV, cD)) = coefs2 #Separando os coeficientes
imgr = pywt.idwt2(coefs2, 'haar', mode = 'periodization')  #1 nível de IDWT

plt.figure(figsize=(10,10))
plt.subplot(2,2,1)
plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")
Out[68]:
Text(0.5, 1.0, 'CD - Bordas Diagonais')

C.1 Cálculo do Erro Quadrático Médio (MSE) e da Relação Sinal Ruído de Pico (PSNR)

a) A MSE é obtida calculando somando-se o erro quadrático de reconstrução pixel a pixel entre a Imagem Original (O) da Reconstruída (R) e normalizando pela dimensão (LxA) da imagem:

$ MSE = \frac{1}{LA}{\sum_{i=0}^L}{\sum_{j=0}^A [O(i,j) - R(i,j)]^2}$

b) A SNR de pico (PSNR) é definida para cada plano componente da imagem como:

$ PSNR = 10.log_{10} \left( \frac{{MAX_I}^2}{MSE} \right) $

sendo $MAX_I$ o valor máximo do pixel, que para 8 bits equivale a 255, logo:

$ PSNR = 20.log_{10}(255) - 10.log_{10} (MSE) $

OBS.: Para uma imagem RGB, $ MSE = MSE_R + MSE_G + MSE_B $, sendo similar definiação para YCrCb e HSV

In [69]:
# Calculo da MSE P&B
A, L, Camadas = img.shape
dif = img_gray - imgr
MSE_gray = np.sum(np.matmul(dif,np.transpose(dif)))/(A*L)
print("MSE_Y = {:.2e}".format(MSE_gray))
PSNR_Y = 20*np.log10(255) - 10*np.log10(MSE_gray)
print("PSNR_Luma = {:.2f} dB".format(PSNR_Y))
plt.figure(figsize=(20,10))
infograf = "Imagem Reconstruída de Luminância (Y) com PSNR = " + str(np.uint8(PSNR_Y)) + ' dB'
plt.subplot(1,2,1); plt.imshow(img_gray,'gray'); plt.title("Imagem Original P&B")
plt.subplot(1,2,2); plt.imshow(imgr,'gray'); plt.title(infograf)

MSE_Y = 9.10e-13
PSNR_Luma = 168.54 dB
Out[69]:
Text(0.5, 1.0, 'Imagem Reconstruída de Luminância (Y) com PSNR = 168 dB')

D. Teste das Funções de Multiresolução wavedec2() e waverec2()

In [70]:
C = pywt.wavedec2(img_gray,'haar', mode = 'symmetric', level=2) # Dois níveis de decomposição DWT
imgr2 = pywt.waverec2(C, 'haar', mode = 'symmetric') # Dois níveis de IDWT

# Para extrair os coeficientes de cada nível
cA2 = C[0]  # Coeficientes de Aprosimação nível 2
(cH1, cV1, cD1) = C[-1] # Coeficientes de Detalhes nível 1
(cH2, cV2, cD2) = C[-2] # Coeficientes de Detalhes nível 2

# Imagem Original
img_gray = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

# Plot dos coeficientes do nível 2
plt.figure(figsize=(5,5))
plt.subplot(2,2,1)
plt.imshow(cA2, 'gray'); plt.title('Ap. N2: CA2')
plt.subplot(2,2,2)
plt.imshow(cV2, 'gray'); plt.title('B. V. N2: CV2')
plt.subplot(2,2,3)
plt.imshow(cH2, 'gray'); plt.title('B. H. N2: CH2')
plt.subplot(2,2,4)
plt.imshow(cD2, 'gray'); plt.title('B. D. N2: CD2')

# Plot Original e Reconstrução
plt.figure(figsize=(20,10))
plt.subplot(1,2,1); plt.imshow(img_gray,'gray'); plt.title('Imagem Original')
plt.subplot(1,2,2); plt.imshow(imgr2,'gray'); plt.title('Imagem Reconstruída')
Out[70]:
Text(0.5, 1.0, 'Imagem Reconstruída')

E. Efetuar uma "Montagem" com wavedec2() e wavedecn()

1º Nível

In [70]:
In [71]:
CV1 = cV1.copy()
CH1 = cH1.copy()
CD1 = cD1.copy()

2º Nível

In [72]:
CA2 = cA2.copy()
CH2 = cH2.copy()
CV2 = cV2.copy()
CD2 = cD2.copy()

# Matriz Final Completa
CC1 = np.bmat([[CA2,CV2],[CH2,CD2]])
h,w = CV1.shape
CC1 = CC1[:h,:w]
CC = np.bmat([[CC1,CV1],[CH1,CD1]])

print(CA1.shape, CH1.shape, CV1.shape, CD1.shape)

(640, 481) (640, 480) (640, 480) (640, 480)
In [73]:
plt.figure(figsize=(20,20))
plt.imshow(CC,'gray')
plt.title('Codificação de Imagem em multinível com função wavedec2()')
Out[73]:
Text(0.5, 1.0, 'Codificação de Imagem em multinível com função wavedec2()')

F. Reconstrução de Imagem Colorida

In [74]:
# Imagem Original

plt.figure(figsize=(20,20))
plt.imshow(img); plt.title("Imagem Original")

# Codificação por planos de cores
# Plano Vermelho
coefs_R = pywt.dwt2(img[:,:,0],'haar', mode='periodization')  #1 nível de DWT R
(cA_R, (cH_R, cV_R, cD_R)) = coefs_R #Separando os coeficientes
cr_R = pywt.idwt2(coefs_R, 'haar', mode = 'periodization')  #1 nível de IDWT R

plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_R,'Reds_r'); plt.title("CA_Red - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_R,'Reds_r'); plt.title("CV_Red - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_R,'Reds_r'); plt.title("CH_Red - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_R,'Reds_r'); plt.title("CD_Red - Bordas Diagonais")

plt.figure(figsize=(20,20))
plt.imshow(cr_R, 'Reds_r'); plt.title("Imagem Reconstruída Red")
Out[74]:
Text(0.5, 1.0, 'Imagem Reconstruída Red')
In [75]:
# Plano Verde
coefs_G = pywt.dwt2(img[:,:,1],'haar', mode='periodization')  #1 nível de DWT G
(cA_G, (cH_G, cV_G, cD_G)) = coefs_G #Separando os coeficientes
cr_G = pywt.idwt2(coefs_G, 'haar', mode = 'periodization')  #1 nível de IDWT G

plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_G,'Greens_r'); plt.title("CA_Green - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_G,'Greens_r'); plt.title("CV_Green - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_G,'Greens_r'); plt.title("CH_Green - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_G,'Greens_r'); plt.title("CD_Green - Bordas Diagonais")

plt.figure(figsize=(20,20))
plt.imshow(cr_G, 'Greens_r'); plt.title("Imagem Reconstruída Green")
Out[75]:
Text(0.5, 1.0, 'Imagem Reconstruída Green')
In [76]:
# Plano Azul
coefs_B = pywt.dwt2(img[:,:,2],'haar', mode='periodization')  #1 nível de DWT B
(cA_B, (cH_B, cV_B, cD_B)) = coefs_B #Separando os coeficientes
cr_B = pywt.idwt2(coefs_B, 'haar', mode = 'periodization')  #1 nível de IDWT B

plt.figure(figsize=(20,20))
plt.subplot(2,2,1)
plt.imshow(cA_B,'Blues_r'); plt.title("CA_Blue - Aproximação")
plt.subplot(2,2,2)
plt.imshow(cV_B,'Blues_r'); plt.title("CV_Blue - Bordas Verticais")
plt.subplot(2,2,3)
plt.imshow(cH_B,'Blues_r'); plt.title("CH_Blue - Bordas Horizontais")
plt.subplot(2,2,4)
plt.imshow(cD_B,'Blues_r'); plt.title("CD_Blue - Bordas Diagonais")

plt.figure(figsize=(20,20))
plt.imshow(cr_B, 'Blues_r'); plt.title("Imagem Reconstruída Blue")
Out[76]:
Text(0.5, 1.0, 'Imagem Reconstruída Blue')
In [77]:
# Reconstrução Nível 1 Colorida 
A1, L1 = cA_R.shape
imgrec1 = np.zeros((A1,L1,3))
imgrec1[:,:,0] = cA_R.copy()
imgrec1[:,:,1] = cA_G.copy()
imgrec1[:,:,2] = cA_B.copy()
plt.figure(figsize=(10,10))
plt.imshow(imgrec1); plt.title("Imagem Reconstruída DWT/IDWT Nível 1")

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[77]:
Text(0.5, 1.0, 'Imagem Reconstruída DWT/IDWT Nível 1')

G. Salvando as Aproximações e depois fazendo download dos arquivos, calcular a taxa de compressão com o original

In [78]:
cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N1_Y.bmp', cA) # Aproximação Nível 1 só Y

cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N2_Y.bmp', cA2) # Aproximação Nível 2 só Y
Out[78]:
True

H. Gravando o Arquivo Codificado DWT/IDWT nível 1 Colorido, calcular a taxa de compressão com o original

In [79]:
cv.imwrite('/content/gdrive/My Drive/UFABC2019/CSM/lab4/bruno_DWT_N1_colorida.bmp', imgrec1) # Gravando Aproximação Nível 1 Colorida
Out[79]:
True

I. Reconstrução da Imagem colorida e Cálculo da MSE de cada plano de cor e da PSNR total

In [80]:
# Reconstrução Colorida Original
A,L = imgr2.shape
imgrec = np.zeros((A,L,3))
imgrec[:,:,0] = cr_R.copy()
imgrec[:,:,1] = cr_G.copy()
imgrec[:,:,2] = cr_B.copy()

# Calculo do MSE colorida
dif2 = img - imgrec
MSE_R = np.sum(np.matmul(dif2[:,:,0],np.transpose(dif2[:,:,0])))/(A*L) # Erro Quadrático Médio plano R
MSE_G = np.sum(np.matmul(dif2[:,:,1],np.transpose(dif2[:,:,1])))/(A*L) # Erro Quadrático Médio plano G
MSE_B = np.sum(np.matmul(dif2[:,:,2],np.transpose(dif2[:,:,2])))/(A*L) # Erro Quadrático Médio plano B
print("MSE_Red= {:.2e}".format(MSE_R), " MSE_Green= {:.2e}".format(MSE_G), " MSE_Blue= {:.2e}".format(MSE_B))

# Cálculo da SNR de pico colorida (PSNR), 3 camadas R, G e B
PSNR = 20*np.log10(255) - 10*np.log10(MSE_R + MSE_G + MSE_B)
print("PSNR total = {:.2f} dB".format(PSNR))

plt.figure(figsize=(20,20))
infograf2 = "Imagem Reconstruída com PSNR total = " + str(np.uint8(PSNR)) + ' dB'
plt.imshow(imgrec); plt.title(infograf2)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

MSE_Red= 1.28e-12  MSE_Green= 8.19e-13  MSE_Blue= 6.65e-13
PSNR total = 163.72 dB
Out[80]:
Text(0.5, 1.0, 'Imagem Reconstruída com PSNR total = 163 dB')

I. Cada um deve repetir todos os passos para a sua foto individual

J. Relatório (página)

  1. Acrescentar também a foto-montagem do grupo todo, mas codificar esta foto-montagem em um nível com DWT também as componentes Cr e Cb (efetuamos aqui apenas para a componente Y)
  2. Fazer uma tabela com as taxas de compressão obtidas e a PSNR de cada uma das fotos utilizadas
In [81]:
# Imagem do grupo
img = cv.imread('/content/gdrive/My Drive/UFABC2019/CSM/lab4/grupo.jpg')

# Convertendo para YCrCb e extraindo canais
img_YCrCb = cv.cvtColor(img, cv.COLOR_BGR2YCrCb)
#Y,Cr,Cb = cv.split(img_YCrCb)
Y,Cr,Cb = img_YCrCb[:,:,0], img_YCrCb[:,:,1], img_YCrCb[:,:,2]
In [82]:
print("Componente Y")
Y_DWT = pywt.dwt2(Y,'haar', mode='periodization')  #1 nível de DWT
(cA, (cH, cV, cD)) = Y_DWT
Y_IDWT = pywt.idwt2(Y_DWT, 'haar', mode = 'periodization')  #1 nível de IDWT

plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")

Componente Y
Out[82]:
Text(0.5, 1.0, 'CD - Bordas Diagonais')
In [83]:
print("Componente Cb")

Cb_DWT = pywt.dwt2(Cb,'haar', mode='periodization')  #1 nível de DWT
(cA, (cH, cV, cD)) = Cb_DWT
Cb_IDWT = pywt.idwt2(Cb_DWT, 'haar', mode = 'periodization')  #1 nível de IDWT

plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")

Componente Cb
Out[83]:
Text(0.5, 1.0, 'CD - Bordas Diagonais')
In [84]:
print("Componente Cr")

Cr_DWT = pywt.dwt2(Cr,'haar', mode='periodization')  #1 nível de DWT
(cA, (cH, cV, cD)) = Cr_DWT
Cr_IDWT = pywt.idwt2(Cr_DWT, 'haar', mode = 'periodization')  #1 nível de IDWT

plt.figure(figsize=(18,8))
plt.subplot(2,2,1);plt.imshow(cA,'gray'); plt.title("CA - Aproximação")
plt.subplot(2,2,2);plt.imshow(cV,'gray'); plt.title("CV - Bordas Verticais")
plt.subplot(2,2,3);plt.imshow(cH,'gray'); plt.title("CH - Bordas Horizontais")
plt.subplot(2,2,4);plt.imshow(cD,'gray'); plt.title("CD - Bordas Diagonais")

Componente Cr
Out[84]:
Text(0.5, 1.0, 'CD - Bordas Diagonais')
In [85]:
# Calculo da MSE P&B
A, L, Camadas = img_YCrCb.shape
imgr = np.zeros(img_YCrCb.shape, dtype=Y_IDWT.dtype)
imgr[:A,:L,0] = Y_IDWT[:A,:L]
imgr[:A,:L,1] = Cr_IDWT[:A,:L]
imgr[:A,:L,2] = Cb_IDWT[:A,:L]

def MSE(A,B):
    difference_array = np.subtract(A, B)
    squared_array = np.square(difference_array)
    mse = squared_array.mean()
    return mse
def PSNR(mse):
    return 20*np.log10(255) - 10*np.log10(mse)

MSE_YCrCb = MSE(img_YCrCb.astype(np.float64), imgr)
print("MSE_YCrCb = {:.2e}".format(MSE_YCrCb))

PSNR_YCrCb = PSNR(MSE_YCrCb)
print("PSNR_YCrCb = {:.2f} dB".format(PSNR_YCrCb))

plt.figure(figsize=(18,8))
plt.subplot(1,2,1); plt.imshow(img_YCrCb,'gray'); plt.title("Imagem YCrCb")
infograf = "Imagem Reconstruída (YCrCb) com PSNR = {:.2f} dB".format(PSNR_YCrCb)
plt.subplot(1,2,2); plt.imshow(np.uint8(imgr),'gray'); plt.title(infograf)

MSE_YCrCb = 1.53e-27
PSNR_YCrCb = 316.29 dB
Out[85]:
Text(0.5, 1.0, 'Imagem Reconstruída (YCrCb) com PSNR = 316.29 dB')